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Monte Carlo methods were used to calculate heat capacities as functions of temperature for 
classical atomic clusters of aggregate sizes 25 < N < 60 that were bound by pairwise Lennard-Jones 
potentials. The parallel tempering method was used to overcome convergence difficulties due to 
quasiergodicity in the solid-liquid phase-change regions. All of the clusters studied had pronounced 
peaks in their heat capacity curves, most of which corresponded to their solid-liquid phase-change 
regions. The heat capacity peak height and location exhibited two general trends as functions of 
cluster size: for N — 25 to 36, the peak temperature slowly increased, while the peak height slowly 
decreased, disappearing by N = 37; for N = 30, a very small secondary peak at very low temperature 
emerged and quickly increased in size and temperature as N increased, becoming the dominant peak 
by N — 36. Superimposed on these general trends were smaller fluctuations in the peak heights that 
corresponded to "magic number" behavior, with local maxima found at N — 36, 39, 43, 46 and 49, 
and the largest peak found at N — 55. These magic numbers were a subset of the magic numbers 
found for other cluster properties, and can be largely understood in terms of the clusters' underlying 
geometries. Further insights into the melting behavior of these clusters were obtained from quench 
i-G ' studies and by examining rms bond length fluctuations. 
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I. INTRODUCTION 

The investigation of small atomic and molecular 
clusters remains an active area of.-j-esearch, both 
theoreticallycTO and experimentallyE-TO Computer 
simulations of clusters have revealed a wide diversity of 
fundamental behavior that is highly sensitive to many 
factors, such as the cluster aggregate size N, the poten- 
tial used to model the intermolecular interactions, and 
even the simulation method itself. The Lennard-Jones 
(LJ) potential is the most commonly used model po- 
tential for computer simulations of cluster properties, 
and it is especially useful for modeling rare gas clusters. 
A key aspect of this model when applied to small and 
medium sized atomic clusters is that the resulting poten- 
tial surfaces have global minima that are dominated by 
icosahedral-likc structures, with fivefold symmetries that 
are markedly different than the closest-packed structure 
found in corresponding bulk systems. 

A major motivation for studying clusters is the in- 
sight that can be obtained in understanding the tran- 
sition from finite to bulk behavior. Because of their 
high proportion of surface atoms, many clusters have 
physical properties that typically exhibit a very irreg- 
ular dependence on their aggregate size, with certain 
sizes standing out in particular. For example, the mass 
spectra of clusters typically exhibit especially abundant 
sizes that often reflect particularly stable structures, .spufe. 
daily reactive clusters, or closed electronic shells Jldliil'tll 
These "magic, number" sizes have been of great theoreti- 
cal interest ,ElrE^I and much work has been done in relat- 
ing magic number sequences to cluster structure in terms 
of "soft" sphere packings, since many magic number sizes 
in LJ clusters such as N = 13, 19, 55, and 147 correspond 
to compact structures that are especially stable.Elu Many 
cluster properties obtained from simulation studies such 
as binding energies, root mean square (rms) bond length 
fluctuations and heat capacities also show very irregu- 
lar dependencies on the cluster size that mostly corre- 
spond to their usual magic number sequences. For exam-, 
pie, the heaLjcaoacity peak values per atom for LJ13&EH 
and LJ5 5 E§L£rLD are especially large, which can be at- 
tributed to the exceptional thermodynamic stability of 
their lowest-energy isomers relative to their next higher 
energy configurations. To more easily compare various 
cluster property size dependencies in this report, I have 
generalized the concept of magic numbers for arbitrary 
cluster properties to refer to those cluster sizes having 
particularly enhanced values relative to their immediate 
neighbors. In a previous study, I examined these gen- 
eralized magic number effects in LJ clusters ranging in 
size from 4 to 24 atomsJl^l In this study, I extend the 
sequence to 60 atoms, which completes the filling of the 
second Mackay icosahedral shell (to TV = 55) and just 
begins the filling of the third shell (to N = 147). 

The issue of cluster phase transitions is fundamentally 
important, and so has been of great interest. Because 



of their small size, clusters do not have the sharp first 
order transition characteristic of bulk melting, but in- 
stead show changes occurring ovgr a range of tempera- 
tures, typically a few degrees KE Although the phase- 
change range generally decreases as the cluster size in- 
creases, coalescing to the bulk transition temperature 
as N approaches infinity, both the range and the lo- 
cation show much variation as a function of N. Clus- 
ter solid-liquid phase-change regions can be determined 
from computer simulations typically by identifying char- 
acteristic changes in certain cluster properties. For ex- 
ample, cluster potential energies often change markedly 
as functions of temperature in the phase-change regions, 
resulting in distinct peaks in the corresponding heat ca- 
pacity curves; rms bond length fluctuations often rise 
abruptly with increasing temperature at the beginning of 
the phase-change regions, and the fractions of quenched 
cluster configurations corresponding to the lowest-energy 
isomers often decrease sharply with increasing tempera- 
ture in these regions. Defining the phase-change region 
on the basis of such behavior is somewhat ambiguous, 
though, since the temperature ranges associated with 
these different properties do not always coincide, and oc- 
casionally do not even overlap. For a given property, 
the phase-change region temperature range can also de- 
pend on the simulation method; for example, for Ari3, 
the rise in the rms bond length fluctuations as calculated 
from canonical Monte Carlo simulations occurred about 
7 K lower than that calculated from molecular dynamics 
simulations .0 Because of these difficulties, it is convenient 
to loosely define the solid-liquid phase-change region as 
the temperature region between the lower-temperature, 
solid-like region (where the clusters are mostly confined 
to their lower energy configurations, undergoing isomer- 
izations only rarely at the higher-temperature end of 
the region) and the higher-temperature, liquid-like re- 
gion (where the clusters exhibit easy isomcrization and 
diffusion). Cluster behavior in the solid- liquid phase- 
change regions is particularly complicated, with some 
cluster sizes showing coexistence behavior, where "solicL 
like" and "liquid- like" isomers dynamically coexist, Em 
while other cluster sizes show a relatively smooth progres- 
sion of isomerizations between locally similar but globally 
distinct configurations occurring throughout this regionrl 
Although clusters having magic number sizes such as 13 
and 55 have pronounced heat capacity peaks that cor- 
respond closely to the solid-liquid phase-change region 
determined from other cluster properties, other cluster 
sizes such as N = 15, 16 and 17 do not have heat ca- 
pacity peaks in their phase-change regions Jld Since heat 
capcity peaks ha,vc. ^tea-been used to identify phase- 
change regions j£j'Ej'E3SoO characterizing their tem- 
perature dependencies for sequences of LJ clusters can 
help elucidate the underlying reasons for such varied be- 
havior. 

The accurate determination of cluster properties such 
as the heat capacity is hampered by the poor convergence 
of standard simulation methods in the phase-change re- 



gions. This poor convergence is a consequence of quasier- 
godici|fcj4, or the incomplete sampling of configuration 
space Jia Various methods have been developed in re- 
cent years to reduce the systematic errors resulting from 
quasiergodicity, including histogram methods,E2l jump- 
walking methods-( J- walking) J2j~E3 smart walking loath-, 
ods (S- walking) ,E2l and parallel tempering methods.E£rE3 
Many of these methods are based on the coupling of 
configurations obtained from ergodic higher-temperature 
simulations to the quasiergodic lower-temperature simu- 
lations. In my previous study of magic number effects in 
cluster heat capacities, I used J- walking to successfully 
overcome quasiergodicity. Monte Carlo J-walking meth- 
ods couple the usual small-scale Metropolis moves made 
by a lower-temperature random walker, with occasional 
large-scale jumps that move the random walker out of the 
confined regions of configuration space that it could not 
otherwise easily escape within the duration of the sim- 
ulation. There are two complimentary implementations 
of J-walking. In tandem J-walking, higher-temperature 
walkers and lower-temperature walkers are run in tan- 
dem, with the higher-temperature walkers occasionally 
feeding configurations to the lower-temperature walkers. 
Because J-walking obeys detailed balance only approxi- 
mately, the walkers become correlated whenever a lower- 
temperature walker accepts a jump, and so it is necessary 
to run several extra Metropolis passes after a jump is ac- 
cepted to break these correlations, thereby decreasing the 
overall efficiency. The second method breaks the correla- 
tions by running the higher-temperature walkers before- 
hand, and storing representative samples of the config- 
urations in external files, which the lower-temperature 
walkers sample periodically in a random fashion. Thus, 
sampling from an external distribution is more efficient, 
provided that the required distributions are not too large, 
and this was the implementation I used in the previous 
study. It was also successful for the first few clusters 
considered in this study, but for clusters ranging in size 
from N = 31 to 37, which had two heat capacity peaks, 
the heat capacity curves obtained with different J-walker 
distributions showed poor reproducibility for the lower- 
temperature peaks. Since the distribution file sizes al- 
ready totalled more than 2 Gb, it seemed prudent to 
abandon that implementation and switch to the tandem 
J-walker impleHientation. However, a recent study by 
Neirotti, et alr3 of LJ38, which is notoriously difficult to 
simulate accurately, showed that the parallel tempering 
method was remarkably successful in overcoming quasier- 
godicity and calculating an accurate heat capacity curve. 
Parallel tempering is nearly identical to a parallel version 
of tandem J-walking, except that instead of copying the 
higher-temperature J-walker's configuration to the lower- 
temperature walker's configuration whenever a jump is 
accepted, the two configurations are exchanged. This has 
the advantage of maintaining detailed balance, thereby 
eliminating the need to run extra Metropolis passes after 
accepting an exchange, which makes the method more 
efficient than tandem J-walking. Thus, I decided to use 



parallel tempering for all of the cluster sizes examined in 
this study. 

I begin in Section NH with a brief description of the com- 
putational methods used in this work and a summary of 
the calculations undertaken. Section III describes the 



results obtained for the cluster potential energies, heat 
capacities, rms bond length fluctuations and quenched 
configuration distributions, as functions of temperature 
and cluster size. In Section |TVJ , I summarize the impor- 
tant conclusions concerning magic number behavior in 
LJ clusters for N < 60. 



II. COMPUTATIONAL METHODS 

The clusters studied were modeled by the pairwise ad- 
ditive Lennard-Jones potential, 
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Because small clusters can dissociate at higher temper- 
atures and make the cluster definition ambiguous, it is 
customary to confine the clusters by a constraining po- 
tential. For classical Monte Carlo simulations, a perfectly 
reflecting spherical constraining potential of radius i? c , 
centered on the cluster's center of mass, is most conve- 
nient, and was used in this work. 

The classical internal energy and heat capacity were 
calculated by the usual expressions for an TV-atom cluster 
(in reduced units), 
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where U* = U/e, V* = V/e, C^ = C v /k B , and T* = 
ksT/e. Since so many rare gas cluster simulations have 
been done for Ar, I found it convenient to generate the 
temperature-dependent curves using Ar parameters, with 
a = 3.405 A and e = 119.4 K. 

Relative rms bond length fluctuations have 
culated in many previous cluster simulations.! 
These were calculated as 
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where the summations are taken over the entire walk. Be- 
cause of the summations being taken over the entire walk, 
rms bond length fluctuations cannot be determined for 
the J-walking and parallel tempering simulations, whose 
configuration jumps and exchanges effectively reduce the 



long Metropolis walks to a collection of very many short 
excursions about the last accepted jump or exchange 
configuration. Thus, quasi-ergodicity remains a problem 
for rms bond length calculations. Caution must also be 
exercised when interpreting S(T) curves, since they are 
dependent on simulation parameters such .a s . th e walk 
length, as well as the ensemble that is usedBE^Jta 



findings imply that J-walking was successful in access- 
ing all of the important regions of configuration space 
for these clusters, even though it was unable to do so in 
the statistically representative manner needed to accu- 
rately reproduce the heat capacities at very low temper- 
atures. Thus, higher temperature J-walker distributions 
were also generated for cluster sizes 39 to 60 to obtain sets 
of the important low-energy quenched configurations. 



A. J-walking 

J-walking runs were done for cluster sizes 25 < N < 38, 
using sampling from externally stored distributions as de- 
scribed in the preceding study.tJ Distributions for each 
temperature consisted of several files, totaling 2.5 x 10 6 
configurations, sampled every 100, 50, 25, or 10 passes, 
depending on the temperature (100 passes for higher 
temperatures having broad distributions, 10 passes for 
very low temperatures having very narrow distributions) . 
The distributions were generated in stages from higher 
temperatures to lower, beginning with temperatures well 
within the cluster liquid region (50 to 60 K). Data accu- 
mulation was done with walks consisting of 2 x 10 6 total 
passes, with jumps attempted randomly with a frequency 
of 10%. Since any systematic errors in the higher tem- 
perature J-walker distributions can feed into the lower 
temperature distributions and corrupt them, heat capac- 
ity peaks occurring at very low temperatures can be very 
difficult to calculate accurately. In previous studies of 
13-atom binary clusters,Efl very low temperature heat ca- 
pacity peaks corresponding to mixing anomalies could 
be obtained accurately by repeating the full J-walking 
temperature scan several times, each time with newly 
generated distributions, and averaging the independent 
curves. This procedure was followed in this work as well, 
and at least five J-walking runs were done for cluster 
sizes 25 to 30, with the heat capacity curves showing very 
good reproducibility in each case. For cluster sizes 31 to 
38 though, the curves showed very poor reproducibility 
for the lower-temperature peaks, and so the method was 
abandoned. 

Despite the J-walking method's inability to reproduce 
the lower-temperature heat capacity peaks accurately, 
the method was still useful for obtaining a listing of the 
important lowest-energy configurations for each cluster 
size. In each case, representative subsets of the distribp-. 
tion files were quenched using steepest descent methodsta 
and all of the unique isomers found were extracted. While 
this crude method cannot identify all of the configura- 
tions for clusters of these sizes, it is well suited for iden- 
tifying most of the lower energy isomers that dominate 
the low temperature behavior of clusters. For example, 
the quenched J-walker distributions for LJ38 contained 
all xtf the lowest-energy configurations found by Doye et 
al.f-i including the global minimum face-centered-cubic 
(fee) truncated octahedron, which is notoriously difficult 
to obtain from randomly initialized simulations. These 



B. Metropolis Monte Carlo 

Metropolis simulations were also run for each cluster. 
These provided a check of the J-walking and parallel 
tempering results for those temperature regions where 
quasiergodicity in the Metropolis runs was not a problem, 
as well as revealing trends in the systematic errors arising 
from quasiergodicity in the Metropolis runs themselves. 
Temperature scans were also generated using Ar param- 
eters, with the temperature mesh size set to AT = 1K, 
and sometimes to 0.5 K in the phase-change regions. For 
each temperature, simulations consisted of 10 5 warmup 
passes, followed by 10 7 passes with data accumulation. 
The scans in each case were started at the lowest temper- 
ature from the global minimum configuration obtained 
from the corresponding J-walker simulations, and were 
continued well past the cluster solid-liquid phase-change 
region, with the final configuration for each temperature 
used as the initial configuration for the subsequent tem- 
perature. To gain a better appreciation of the extent of 
quasiergodicity in the solid-liquid phase-change region, 
additional Metropolis simulations consisting of 10 8 to- 
tal passes per temperature point were run for each clus- 
ter over a narrower temperature range encompassing the 
phase-change region. 



C. Parallel Tempering 

Parallel temperingj23~E3 is very closely related to 
the tan dem . and parallel walker implementations of J- 
walking,Oc3 with the essential difference being that 
instead of the higher-temperature configuration being 
copied over the lower-temperature configuration when- 
ever a jump is accepted, the two configurations are ex- 
changed. This close similarity between parallel temper- 
ing and tandem J-walking makes it very easy to con- 
vert a tandem J-walking simulation program to a paral- 
lel tempering one. The parallel tempering method's con- 
figuration exchange ensures that it obeys detailed bal- 
ance (provided that the simulation has been warmed 
up enough to be in the asymptotic region), so there 
is no need to run the extra Metropolis passes to break 
correlations between the higher-temperature and lower- 
temperature walkers whenever an exchange move is ac- 
cepted, which is required for tandem J-walking. In ad- 
dition to offering greater efficiency, the obeying of de- 



tailed balance allows parallel tempering simulations to 
vary the exchange-attempt frequency in temperature re- 
gions where the exchange-acceptance rate is too low. In 
J-walking simulations, the jump-attempt frequency was 
typically kept at 10%, since larger values required even 
more extra Metropolis passes to be run whenever a jump 
was accepted to break the correlations, thereby offsetting 
the efficiency gains made by increasing the jump-attempt 
frequency in the first place. 

As with J-walking, the number of temperatures used 
in a parallel tempering simulation and their spacing can- 
not be chosen arbitrarily, but must be selected to en- 
sure that the exchanges are accepted sufficiently often. 
Generally, more temperature walkers and smaller spac- 
ing are required for the solid-liquid phase-change region, 
and for very low temperatures. Since the heat capacity 
curves as a functions of temperature typically had peaks 
in the solid-liquid phase-change regions whose location 
and height I wanted to determine accurately, the num- 
ber of temperatures and their spacing was also affected by 
the need to properly represent the heat capacity curves 
with a sufficient density of points to allow for good inter- 
polation. This latter requirement resulted in my choos- 
ing temperature spacings in the phase-change region that 
were small enough that the exchange acceptance ratios 
were typically quite high. For very low temperatures, the 
classical heat capacity curve is nearly linear, thus requir- 
ing fewer points to represent. However, the exchange ac- 
ceptance ratio decreases rapidly at lower temperatures, 
requiring more temperature walkers. Since the overall 
computation time is nearly proportional to the number 
of walkers, I was able to reduce the number of lower- 
temperature walkers needed by increasing the exchange- 
attempt frequency Px at lower temperatures, according 
to the relation 
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where P X (T L ) = (1 + 3P x (T H ))/4. This expression 
increases the lowest-temperature exchange-attempt fre- 
quency Px(Tl) by a quarter of the range above the base- 
line highest-temperature exchange-attempt frequency, 
Px(T H ), which was set to 10%. 

The parallel tempering simulations done in this study, 
were similar to those done by Neirotti et al. for LJ38.E3 
For each cluster simulation, a number of random walkers 
at various temperatures were run in parallel, with each 
temperature walker occasionally attempting an exchange 
move with the next higher temperature walker, except 
for the highest temperature walker, which attempted no 
exchanges directly (but did exchange configurations ac- 
cording to attempts made by the next-highest tempera- 
ture walker). Since I had the lowest-energy isomers for 
each cluster from the J-walking distribution quenches, 
for most of the clusters I initialized the temperature 
walkers by using the lowest-energy isomer for the lowest- 
temperature initial configuration, the next lowest-energy 



isomer for the next lowest-temperature initial configura- 
tion, and so on. For the other clusters, I just used the 
lowest-energy isomer to initialize all of the initial config- 
urations. It did not seem to matter which initialization 
scheme was used, and of course, random configurations 
could have been used, but at the cost of longer warmup 
periods. Because I initialized the simulations from the 
lowest-energy isomers, I could quickly run preliminary 
simulations that gave qualitatively good results, which I 
could use to fine-tune the number of parallel tempering 
temperature walkers and their spacing. These prelimi- 
nary parallel tempering runs consisted of 10 5 Metropolis 
warmup passes for each temperature walker, followed by 
10 5 parallel tempering warmup passes and then by 10 
parallel tempering walks with preliminary data accumu- 
lation, with each walk being 10 5 passes in length. The 
configurations for each temperature walker were saved at 
the end of the warmup period, and at the end of each 
walk. The heat capacity curves for each walk were ex- 
amined for consistency, and if significant variations were 
found, the parallel tempering run was repeated (without 
the Metropolis warmup passes) from the last saved con- 
figurations, so that in effect, all of the preceding parallel 
tempering runs served as warmup runs for the final run. 
For some of the clusters, only parallel tempering warmup 
passes were run, indicating that the preceding Metropolis 
warmup passes were likely not necessary (although their 
computational expense was minor). After the prelimi- 
nary runs appeared to be sufficiently warmed up, a longer 
parallel tempering run of 10° warmup passes followed by 
10 7 total passes of data accumulation per temperature 
point was run (10 walks of 10 6 passes each). Again, the 
heat capacity curves for each walk were checked for con- 
sistency; for a few clusters, the variations in the curves 
were still too large, and so the run was repeated from 
the last saved configurations. For LJ38, an especially 
long warmup was, required, consistent with the findings 
of Neirotti et alE3 

One important difference between my implementation 
of parallel tempering and that of Neirotti et al concerns 
the choice of the constraining radius R c . Ideally, the 
value of the constraining radius should be large enough 
that the cluster properties being simulated are not ad- 
versely affected by the constraining potential, but not so 
large that cluster dissociations lead to overly long conver- 
gence times. In my previous study of clusters of size 4 to 
24, I used R c — 4.0 a, a value large enough not to affect 
the heat capacity curves, except at higher temperatures, 
well past the solid-liquid phase-change regions. Thus, 
for this study, I also set R c = 4.0 c for the Metropolis 
and J- walker simulations. Neirotti et al used a value of 
R c = 2.25(7 for their simulation of LJ38, since they had 
difficulties attaining ergodicity with larger constraining 
radii. For LJ38, such a small constraining radius was not 
a problem since the heat capacity peak occurred at a rela- 
tively low temperature. However, for most of the cluster 
sizes included in this study, such a small value would 
affect the heat capacity peak height and location. After 



some preliminary investigations, I chose to set R c = 3.5 a 
for parallel tempering simulations of cluster sizes 25 to 48, 
and to R c — 4.0 a for the larger cluster sizes. Figure [l] il- 
lustrates the effects of the constraining radius on clusters 
Ar3o and Ai'32. The heat capacity peak for R c = 3.0 cr 
is slightly truncated, but the peak for R c = 3.5 a is in 
good agreement with the R c — 4.0 c peak, even though 
the high-temperature side of the peak is substantially 
affected. These two examples also show how well paral- 
lel tempering was able to overcome the quasiergodicity 
evident in the Metropolis curves. The two parallel tem- 
pering curves for Ar32 show very good agreement for the 
smaller peak, which is completely absent in the Metropo- 
lis simulations. 

The converged heat capacity curve I obtained for Ar38 
was mostly similar to the curve obtained by-Neirotti et al. 
in their parallel tempering study of LJ38,t3 with a small 
shoulder in the low-temperature side of the peak. The 
reduced peak temperature 0.1653 was also in good agree- 
ment with that found by Neirotti et al. (0.166), although 
my peak value was slightly higher ({Cv)/NkB = 5.858 
compared to 5.62) and the higher-temperature part of 
my curve rose substantially higher; these discrepancies 
are mostly due to Neirotti et al. having used the smaller 
constraining radius. As pointed out by Neirotti et al., 
the shoulder in the heat capacity curve is in contrast to 
the small lower-temperature peak-found by Doyle, Wales 
and Miller in their study of LJ3§E3 (their larger peak lo- 
cation and height were similar to the parallel tempering 
results, though, at 0.17 and 5.68, respectively). 

Finally, follow-up parallel tempering simulations of 10 6 
passes per temperature were run for each cluster with pe- 
riodic quenching of the configurations to determine the 
population distributions of the lower-energy isomers as 
functions of temperature. These were initialized using 
the final configurations from their corresponding preced- 
ing parallel tempering runs. Since I was primarily inter- 
ested in using the quench results to provide qualitative 
insight into the heat capacity results, rather than deter- 
mining the quench curves themselves with high accuracy, 
these shorter walk lengths were sufficient. 



III. RESULTS 

A. Structural properties 

Much of the magic number behavior of cluster heat ca- 
pacities can be understood in terms of the magic number 
behavior of the underlying potential surfaces. For exam- 
ple, for the small clusters I studied earlier ,Eil the largest 
heat capacity peaks were observed for sizes N — 13, 19 
and 23, corresponding to single, double and triple icosa- 
hedra, respectively. The same magic numbers were found 
in some corresponding structural properties as well, such 
as the binding energy differences, 



AE b (N) = -[V min (N) - V min (N - 1)]. 



(6) 



The lowest-energy configurations for each of the clusters 
I examined in this study were obtained from quenches of 
J-walker distributions. These were consistent with .ike. 
lowest-energy configurations found in other studies JjuEj 
From N = 25 to 30, the global minima all consisted of 
a N = 19 double icosahedral core with atoms added to 
tetrahedrally bonded face sites (anti-Mackay overlayers), 
while for N = 31 to 54 (except for N = 38), the added 
atoms occupied sites corresponding to the outer shell 
the N = 55 Mackay icosahedron (Mackay overlayers); 
N = 38 did not Jit-thc pattern, it being a truncated octa- 
hedron instead. e3H Clusters for N = 56 to 60 had added 
atoms that contributed to the anti-Mackay overlayer of 
the 55-atom core, beginning the building up of the next 
Mackay icosahedron (N = 147). 

Although the lowest-energy isomers play a significant 
role in cluster thermodynamics, other low-energy iso- 
mers also play important roles, and the size of the gap 
between the lowest-energy isomer and the next lowest- 
energy isomer is an esuecially important predictor of 
magic number behavior.aQu Figure @ shows the energy 
levels of the lower-energy configurations for each clus- 
ter size, relative to their global minimum. Sizes N = 
26, 29, 32, 36, 39, 43, 46, 49, 55 and 58 have relatively large 
gaps between their global minimum and next higher min- 
imum, compared to their immediate neighbors, and so 
can be considered a magic number sequence; the gaps 
between the lowest-energy and next lowest-energy iso- 
mers are also plotted in Figure as functions of N. The 
Mackay icosahedron N = 55 has an especially large gap. 
There is also much variation in the densities of the low- 
lying isomers evident in Figure |2j, with densities being 
substantially greater for sizes N — 28 to 37 than for the 
other sizes, while sizes N = 52 to 56 have particularly 
low densities. 

The binding energy differences for the clusters are also 
displayed in Figure |3j. Except for N — 38, the max- 
ima in the binding energy differences are the same as 
the maxima in the energy differences between the two 
lowest-energy isomers. Also shown in Figure H is a plot 
of another indicator of structural magic number behav- 
ior, the second finite difference of the energy, 

A 2 E{N) = V min {N + 1) + V min (N - 1) - 2V min (N). (7) 

These magic numbers are also in agreement with the 
magic numbers obtained for the binding energy differ- 
ences, again, except for N = 38. The discrepancy in the 
magic number sequence for N = 38 is due to that clus- 
ter being unique in this range by having an fee structure 
as its global minimum, so that those measures that are 
relative to the cluster's immediate neighbors are consid- 
erably different than those of the other clusters, whose 
immediate neighbors are structurally similar. 

Figure shows a generalization of the zero- 
temperature A%E(N) curve shown in Figure |3|, with 
the curves extended to non-zero temperatures to form 
a A%E(N,T) surface. The curves were obtained from 



parallel tempering simulations; the temperature scale is 
for Ar. The surface can be seen to encompass two general 
regions: a flat region at higher temperatures, correspond- 
ing to the clusters' liquid-like regions, and a corrugated 
region at lower temperatures, corresponding to the clus- 
ters' solid-like regions. The plot shows that many of the 
clusters are locked in their lowest-energy isomers for ex- 
tended temperature ranges, especially the magic number 
sizes N = 26, 29, 39, 43, 46, 49 and 55 — the A 2 E peak 
values differ very little from their zero-temperature values 
well into the phase-change region, until their abatement 
over a relatively narrow temperature range. The plot also 
shows that the phase-change temperature regions have a 
strong dependence on the cluster size, with the melting 
temperature generally decreasing from about 30 K for 
N = 26 to about 20 K for N = 36, and then generally 
increasing to a maximum of about 40 K for N = 55. 



B. Thermodynamic properties 

Figure H plots the reduced potential energy per particle 
for each cluster size over a temperature range spanning 
the solid-like and liquid-like regions. Magic number be- 
havior can be seen for sizes N = 26, 29, 32, 36, 39, 43, 
46, 49 and 55 by the relatively large differences at lower 
temperatures between their potential energies and those 
of their preceding cluster size. The Mackay icosahedron 
N = 55 is so stable relative to its neighbors that its 
curve actually crosses several of the other curves. The 
phase regions are apparent in Figure by the transition 
from the irregular spacing between the curves seen for 
the lower temperature, solid-like regions, to the regular 
spacing seen for the higher temperature, liquid-like re- 
gions. As was seen in the A2E(N, T) surface in Figure 0, 
the melting temperature generally decreases as the clus- 
ter size increases from N = 25, reaching a minimum of 
about 20 K for N about 36, and then increases to a max- 
imum for TV = 55. 

Curves of the reduced heat capacity per particle for 
each cluster size are shown in Figures |3[ and |S|, while 
Figure shows the corresponding peak height and tem- 
perature as functions of the cluster size; the peak pa- 
rameters are also listed in Table D. The peak parameters 
were obtained by smoothing and interpolating the com- 
bined parallel tempering data (10 walks of 10 6 passes 
per temperature point each). The associated parame- 
ter uncertainties were estimated by determining the peak 
height and location for each of the f walks separately, 
and calculating the standard error of the set of these val- 
ues. As a check, additional runs of f 7 total passes per 
temperature were done for a few of the clusters, and the 
standard deviations of the resulting set of peak param- 
eters were compared to the corresponding uncertainty 
estimates. These were found to be in good agreement. 
The extra check runs also showed the parallel tempering 
heat capacity results to be very reproducible, unlike the 



corresponding J- walking results. 

Each of the clusters had at least one peak in its heat 
capacity curve, and overall, while their was much varia- 
tion in the curves, systematic trends are evident. Clus- 
ter sizes N — 36,39,43,46,49 and 55 had peak heights 
that were substantially greater than those of their imme- 
diate neighbors, and thus represent the magic numbers 
for the heat capacity sequence. These cluster sizes were 
also magic numbers for the A2-E(iV) sequence and the 
lowest-energy isomer gaps shown in Figure pi but the 
other structural magic numbers (N = 26, 29 and 32) did 
not show magic number behavior in their heat capacity 
curves. Instead, for cluster sizes 25 to 36, both the heat 
capacity peak height and the peak temperature showed 
remarkably smooth variation, with the peak temperature 
slowly increasing, while the peak height slowly decreased 
in a nearly linear fashion, until the peak disappeared at 
N = 37. As can be seen in Figure pi this behavior was 
somewhat affected by the choice of the constraining ra- 
dius. For the Metropolis simulations run at R c — 4.0 a, 
the heat capacity curves had greater slopes at higher tem- 
peratures, resulting in the peak disappearing by N = 34. 

More interesting was the emergence and evolution of 
a second heat capacity peak, as the first peak gradually 
disappeared. Beginning with N = 30, a small, lower- 
temperature peak can be seen in Figure |6| to develop 
and gradually increase in both height and temperature, 
until it becomes the dominant peak at N = 36. This 
trend can be seen to continue in Figures p] and @, until 
the peak reaches its greatest height and highest temper- 
ature for the magic number size N — 55. This behavior 
suggests partitioning the heat capacity peaks into two se- 
quences, with the first sequence containing the small-iV 
peaks that evolve into the large, broad high-temperature 
peaks as the second Mackay shell initially builds up, 
while the second sequence originates from the smaller, 
narrower, lower-temperature peaks that evolve into the 
dominant peaks as the second Mackay shell is completed. 
The possible beginning of a third sequence correspond- 
ing to the initial building up of the third Mackay shell 
can be inferred from the curves in Figure g, where Arss 
has a small, lower-temperature peak, and Ar 5 g and Ar 60 
have discernable shoulders in the low-temperature sides 
of their peaks. 

The emergence of the smaller peak nearly coincides 
with the transition at N = 31 from the anti-Mackay over- 
layer in the lowest-energy configuration to the Mackay 
ovcrlayer, and the cluster sizes where two prominent heat 
capacity peaks are evident correspond to those sizes hav- 
ing the greatest density of low- lying isomers (Figure y). 
Doye et al. found that for clusters in this size range, the 
high density of isomers was due to the small differences 
in energy between the various configuration types, anti- 
Mackay, Mackay, decahedral, fee, etc;E3 similar findings 
were obtained for clusters bound by a Morse potential 
whose range parameter is similar to the Lennard- Jones 
potential. El 

The two heat capacity sequences make the determi- 



nation of the cluster melting temperature from heat ca- 
pacity peak temperatures ambiguous for sizes N = 31 
to 36. Peaks in heat capacity curves have typically 
been associated with qlu a tqr phase changes, particularly 
solid-liquid changes. E3E3l23'EaE3Ej However, peaks can 
arise fro m o th er circumstances as well, such as mixing 
anomaliesEJEJ and premelting effects, which can make 
the assignment of the cluster melting temperatures diffi- 
cult. Calvo and Spiegelmann defined the melting temper- 
atures of the sodium clusters they studied whose heat ca- 
pacity curves had multiple peaks,, as the temperature as- 
sociated with the highest peak.123 Applying such a defini- 
tion to the clusters shown in Figure @ that have two heat 
capacity peaks would imply that the reduced melting 
temperature would jump from T* — 0.3784 for N = 35 
to 0.1427 for N = 36. Given the structural similari- 
ties of these two clusters, such a large difference in the 
melting temperature would be surprising, and as will be 
seen later when other cluster properties are examined, 
the two clusters do in fact have similar melting temper- 
atures. Also, the potential energy curves shown in Fig- 
ure and the A 2 E(N, T) surface shown in Figure both 
suggest that the melting temperature for these clusters is 
about T* — 0.17, and thus intermediate to the two heat 
capacity peak temperatures. 

There is much evidence of systematic error due to 
quasiergodicity in the Metropolis results, especially for 
sizes N = 31 to 38, where the Metropolis heat capac- 
ity curves had lower-temperature peaks that were either 
missing (N — 31 to 34) or poorly formed (N = 35 to 
38). For the other cluster sizes, the Metropolis heat 
capacity curves showed very good agreement with the 
parallel tempering results, except for temperatures near 
the heat capacity peak, where the Metropolis results for 
10 7 total passes per temperature point were substan- 
tially lower than the parallel tempering results, and even 
the Metropolis results for 10 8 total passes were typically 
slightly lower. The fact that parallel tempering sim- 
ulations of 10 7 passes could provide significantly bet- 
ter results than similar Metropolis simulations of 10 8 
passes, for not much more computational overhead than 
Metropolis simulations of 10 7 passes, clearly shows the 
superiority of parallel tempering over the Metropolis 
method. 

There have been several previous-. studies of 55-atom 
cluster heat capacities done.tSES&LJLZl The study most 
directly comparable to the parallel tempering results I 
obtained for N — 55 is a J-walking study of Arss done by 
L6pez.c3 He obtained peak parameters of (Cy) /Nks — 
19.5 and T* = 0.29, compared to my parallel tempering 
results of 17.45 and 0.294, respectively. Lopez used the 
same constraining volume radius as I did for N = 55, 
R c = 4.0cr, but his J-walker distributions consisted of 
only 4 x 10 configurations, which are considerably fewer 
than the 2.5 x 10 6 configurations I used for my J-walking 
studies of cluster sizes N — 25 to 38. Given the poor re- 
sults I obtained with J-walking using larger distributions 
for smaller cluster sizes, I suspect that the discrepancy 



in the heat capacity peak parameters is due to Lopez 
having used insufficiently representative J-walker distri- 
butions in his study. The heat capacity curve obtained 
from a Metropolis simulation of 10 8 total passes per tem- 
perature that is shown in Figure @ had peak parameters 
of 17.00 and 0.294, respectively, which are close to the 
parallel tempering results. The small discrepancy be- 
tween these Metropolis results and the parallel temper- 
ing results is similar to the discrepancies obtained for the 
other clusters of similar size. Labastie and Whetten ob- 
tained peak parameters of 15.4 and 0.3 for Ai'55 using 
histogram methods, but their smaller peak height was 
most likely due to their constraining radius being consid- 
erably smaller at R c = 2.6 a. Doye and Wales obtained 
peak parameters of 12.4- and 0.298 for Ai'55 using the su- 
perposition method,! 3 3 H 3 1 which is in very good agreement 
with my results. 



C. Bond length fluctuations 

The rms bond length fluctuation 8 is another measure 
that has been used often in cluster—simulations to esti- 
mate the melting temperature. BqOEjo According to 
the Lindemann criterion, melting occurs in bulk matter 
when 8 is about 10% or greater. However, much care 
must be exercised in interpreting rms bond length fluctu- 
ations obtained from cluster simulations. First, as Calvo 
and Spiegelmann have pointed out, the Lindemann cri- 
terion was originally proposed for bulk systems, imply- 
ing that it does not necessarily extend to finite clusters, 
and the choice of 10% for the critical value is somewhat 
arbitraryO Second, 8{T) is dependent on the simulation 
walk length, with the sharp rise in the curve generally 
occurring |-a.t lower temperatures as the walk length is 
increased." For cluster simulations of very long length, 
a value of 8 w 0.1 can be obtained at temperatures where 
the cluster spends most of its time in some solid-like form, 
with only occasional isomerizations to other solid-like 
forms, rather than undergoing the ready diffusion typical 
of liquid-like behavior. Because 8 is obtained from aver- 
ages over entire Metropolis walks or molecular dynam- 
ics trajectories, it cannot be estimated from simulations 
where configuration scrambling occurs, such as J-walking 
and parallel tempering, and so quasiergodicity issues fur- 
ther complicate matters. Using very long walk-lengths to 
reduce systematic errors arising from quasiergodicity can 
result in S(T) curves whose transitions occur well before 
the melting phase-change. 

Despite these caveats, rms bond length fluctuation 
curves have qualitative validity and can provide useful in- 
sights, especially in conjunction with other cluster prop- 
erties. Curves of 5(T) for each cluster size are shown in 
Figures [l(j, [Ly and [12] for averaged walk-lengths of 10 6 
and 10 7 passes per temperature, while Figure O shows 
the TV-dependence of the temperatures corresponding to 
the Lindcmann-likc threshold 8 = 0.2. The threshold 



temperatures were estimated by smoothing and interpo- 
lating the 8(T) curves, and then searching for the inter- 
polated temperature corresponding to S = 0.2. I chose 
to set the threshold slightly higher than the usual Lindc- 
mann criterion of S — 0.1, since some of the S(T) curves 
had too much noise near S = 0.1 to obtain a reliable esti- 
mate of the temperature corresponding to that threshold, 
and because the threshold is somewhat arbitrary anyway; 
the value 5 = 0.2 was also typically near the inflection 
point of the smoothed curves' transition regions. 

The S(T) curves shown in Figures [HI to O are mostly 
similar, but there are some interesting differences. Gen- 
erally, the steepness of the rise in 6 in the transition 
region increases with cluster size, with the magic num- 
ber sizes TV = 43,46 and 55 having especially sharp 
rises. Evidence of quasiergodicity can be clearly seen 
in some curves: the S(T) curves for N = 31 each have a 
noisy shoulder just before the transition region, while the 
curves for N = 32 each have a small peak; the curves for 
N = 38 show the cluster to be locked in its global mini- 
mum for all of the lower-temperature Metropolis simula- 
tions of length 10 7 passes or fewer, until finally escaping 
at T = 19 K. Evidence of pre-melting in the anti-Mackay 
overlayer of the third Mackay shell can be seen in sizes 
N = 56 to 60 by the emergence and evolution of a shoul- 
der in the lower-temperature part of the S(T) curves. 

Also indicated on the temperature axes in Figures M 
to 12 are the corresponding heat capacity peak temper- 



atures obtained from the parallel tempering simulations 
(likewise, Figures |6| to |S| have the corresponding temper- 
atures associated with the Lindemann-like threshold of 
S = 0.2 indicated on their axes). These clearly show that 
the heat capacity peaks do not necessarily correspond 
to a cluster phase change. For cluster sizes N — 25 to 
29, the heat capacity peak temperatures are well within 
the liquid-like region, according to the S(T) curves, while 
for cluster sizes N = 30 to 36, where two heat capacity 
peaks are evident, the higher-temperature heat capacity 
peak is also clearly in the liquid-like region. Moreover, 
these temperatures are high enough that quasiergodicity 
in these regions of the S(T) curves is not an issue. Like- 
wise, most of the lower-temperature heat capacity peaks 
are found to be in the solid-like region of the 6(T) curves, 
well before the transition region (although ergodicity is- 
sues regarding the 5(T) curves in these regions could be 
involved). Only for sizes N = 35 and TV = 36, are the 
lower-temperature heat capacity peak temperatures simi- 
lar to the S(T) transition regions. For the remaining clus- 
ter sizes shown in Figures [Tfl and O, however, the heat 
capacity peak temperatures correspond quite closely to 
the end of the S(T) transition regions, where the curves 
quickly change slope to mark the start of the liquid-like 
region. 

Magic number behavior in the rms bond length fluctua- 
tions is evident in Figure O, especially for the threshold 
temperature values obtained from the longer Metropo- 
lis simulations of 10 7 passes. The threshold tempera- 
tures are significantly higher for the magic number sizes 



N = 26,29,32,36,39,43,46,49 and 55, relative to their 
immediate neighbors. These are the same magic numbers 
seen in Figure [| for A2-E(7V) and for the energy gaps be- 
tween the lowest-energy and next lowest-energy isomers, 
except for N — 58. The 5 data also account for the ab- 
sence of magic number behavior in the heat capacities for 
sizes N = 26, 29 and 32. Since the large, broad, higher- 
temperature heat capacity peaks seen slowly diminish- 
ing in Figures || and ^ are located well past the cluster 
solid-liquid phase-change regions, whatever magic num- 
ber effects might exist for N — 26, 29 and 32 are masked 
by the systematic evolution of the heat capacity peak it- 
self with increasing cluster size. Classical cluster heat 
capacities are inherently a measure of the width of the 
potential energy distributions. Typically, for tempera- 
tures near the cluster melting temperature, the distribu- 
tions are much wider than those at lower temperatures 
(which span mostly lower-energy, solid-like forms) and 
those at higher temperatures (which span mostly higher- 
energy, liquid- like forms), since these distributions span 
both types of forms. However, the distribution width 
in the liquid region also typically increases with tem- 
perature, since the number of accessible configurations 
also increases. Thus, depending on the energy distribu- 
tion of the higher-energy isomer potential minima, the 
inverse dependence of the heat capacity on temperature 
in Eq. (pj) could result in a peak occurring in the liquid- 
like region of the curve that would have little relation to 
the cluster's solid-liquid phase change. 

The threshold temperatures corresponding to S = 0.2 
in Figure [j~3 also imply that the cluster melting tem- 
perature generally decreases as the cluster size increases 
from N — 25, reaching a minimum near N — 37, then 
generally increases with cluster size, reaching a maxi- 
mum for N = 55, and decreases thereafter; magic num- 
ber effects superimposed on these general trends account 
for the local variations. These general trends are con- 
sistent with those inferred from the A2-E(iV, T) surface 
in Figure |4| and the potential energy curves in Figure |5| 
Casero and Soler calculated melting temperatures for LJ 
clusters of size N = 4 to 34 in the microcannonical en- 
semble, defining the melting temperature to be the tem- 
perature where the free-,energies of the solid and liquid 
clusters become equal.E3 Their results indicate a gener- 
ally increasing trend with cluster size, with reduced melt- 
ing temperatures for sizes N = 25 to 34 all above 0.3. 
These results are very similar to the temperatures seen 
for higher-temperature heat capacity peak temperatures 
in Figure 0, and thus are contrary to the melting temper- 
atures implied by the rms bond length fluctuation data 
shown in Figure O. Since the melting temperatures de- 
duced from the S data are in good agreement with the 
melting temperatures obtained from the heat capacity 
peak temperatures for cluster sizes greater than about 
N = 35, and these temperatures are well below 0.3 for 
N < 40, it is likely that the reduced melting tempera- 
tures for sizes N = 25 to 34 are also well below 0.3. 



D. Quench results 

The temperature dependence on the distribution of 
quenched isomers is shown in Figure [L4| for those Ar clus- 
ters whose heat capacities showed magic number behav- 
ior, and in Figures [15] to [17] for some selected non-magic 
number Ar clusters. Each plot shows the percentage of 
quenches to each of the three lowest-energy isomers, as 
well as the total percentage of quenches to the remaining 
higher-energy isomers; the data for all of these plots were 
all obtained from parallel tempering simulations. Much 
of the behavior exhibited by these curves can be ratio- 
nalized in terms of the energy spacings and densities of 
the lower-energy isomers shown in Figure 0. 

The plots for the magic number sizes are all very sim- 
ilar in form, with quenches in the lower-temperature, 
solid-like region represented solely by the lowest-energy 
isomer, and quenches to the other isomers only appear- 
ing at temperatures near the phase-change region, and 
higher. This is a consequence of the large gap between 
the lowest-energy and next lowest-energy isomers that is 
typical of magic number clusters. The only significant 
differences between the different cluster sizes is the tem- 
perature where the fraction of quenches to the lowest- 
energy isomer drops off substantially, which increases 
with cluster size, and can be seen to be very similar to 
the heat capacity peak temperature. 

The plots for some non-magic cluster sizes shown in 
Figure Il5] for TV < 35 are considerably different. The 
curves for TV — 25 have a similar form to the magic 
number curves in Figure |lj, but both the decrease in 
the fraction of quenches to the lowest-energy isomer and 
the increase in the fraction of quenches to the higher- 
energy isomers are much more gradual. This gradual 
transition is consistent with the more gradual transition 
for the S(T) curve, seen in Figure pj. Also, the tran- 
sition region occurs well below the heat capacity peak 
temperature (but near the transition region indicated by 
S(T)). Ar27 has two very closely spaced lowest-energy 
isomers that are well below its next-higher energy iso- 
mers, and most of the lower-temperature quenches arc to 
the second lowest-energy isomer. The temperature where 
the fraction of quenches to the second lowest-energy iso- 
mer exceeds the fraction to the lowest-energy isomer is 
similar to the peak temperature of the very small lower- 
temperature heat capacity peak. Again, the tempera- 
ture region where the fraction of quenches to the higher- 
energy isomers increases substantially corresponds well 
to the transition region in the 8(T) curves, and is well 
below the higher-temperature heat capacity peak tem- 
perature. The quench curves for TV — 31 and 32 are sim- 
ilar to the TV — 27 curve, with the fraction of quenches 
to the lowest-energy isomer dropping rapidly at very low 
temperatures corresponding closely to the small heat ca- 
pacity peak temperatures, and the fraction of quenches 
to higher-energy isomers gradually increasing at temper- 
atures well below the peak temperature of the larger, 



higher-temperature heat capacity peaks. For Ar3i, the 
density of lower-energy isomers is so great that quenches 
to higher-energy isomers dominate at temperatures even 
below the transition region for the S(T) curves shown in 
Figure nfl. The lower-temperature shoulders appearing 
in the TV = 31 S(T) curves and the behavior of the corre- 
sponding quench curves are indicative of pre-melting and 
an extended phase-change region. Ar32 is a magic num- 
ber according to such structural properties as AEb(N) 
and A 2 E(N), but not according to its heat capacity peak 
parameters. Its quench curves are also different than the 
typical magic-number quench curves seen in Figure \L4. 
The TV = 32 6(T) curves show unusual peaks in the solid 
region, just before their sharp rise. This region just fol- 
lows the lower-temperature heat capacity peak, and ac- 
cording to the quench results, almost all of the quenches 
are to the second lowest-energy isomer. This indicates 
that the peaks in the S(T) curves are a result of quasier- 
godicity (the Metropolis heat capacity curves in Figure g 
also show small, spurious peaks in this region). The 
quench curves for TV = 33 and 34 indicate easy access 
to the other low-lying isomers at low temperatures near 
the small heat capacity peak temperatures. 

Figure lty shows quench results for some non-magic 
number clusters for TV = 37 to 48. Ar37, like Ai'31 and 
Ar34, has an especially high density of low- lying isomers, 
and its quench curves also show a substantial fraction of 
quenches to these isomers at temperatures corresponding 
to the heat capacity peak temperature. For TV = 38, 
the temperature region where the fraction of quenches to 
the lowest-energy fee isomer drops off and the fraction 
of quenches to the lowest-energy icosahedral-like isomer 
(the second lowest-energy isomer) increases corresponds 
to the region just before the heat-capacity peak where a 
shoulder can be seen in Figure ffl, which is consistent with 
other studies of this cluster £3 The quench curves for TV = 
41, 45 and 48 are similar to those for TV = 37, since these 
clusters all have accessible lower-energy isomers that are 
only slightly higher than the lowest-energy isomer, but 
the curves for TV = 42 are similar to the curves for the 
magic number clusters in Figure |l4|, indicating that its 
low-lying isomers are not very accessible. 

The isomer energy levels for TV = 50 are similar to 
those for TV = 45 and 48, and the corresponding quench 
curves shown in Figure [Tfl are likewise similar to those for 
TV = 45 and 48. Ars2 and Ars3 have some isomers hav- 
ing closely spaced potential energies only slightly above 
the lowest-energy isomer, and then very few isomers un- 
til AVinin ~ 2.5. This type of isomer energy distribu- 
tion is similar to the distributions seen for binary Ar-Kr 
clusters, where many mixed isomer configur.a±ions had 
energies similar to the lowest-energy isomer,Ej and the 
quench results are similar too, with the low-lying iso- 
mer quench curves showing plateaus throughout much of 
the solid-like temperature region. However, neither Ar52 
nor Ar53 exhibited the small, very low temperature heat 
capacity peaks that were evident in the binary cluster 
curves. Ar56 has two low-energy isomers corresponding 
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to the two unique ways a lone atom can be positioned 
on a 55-atom Mackay icosahedral core, separated by a 
large energy gap from the other isomers corresponding 
to rearrangements of the 55-atom core. The TV = 56 
quench curves show substantial fractions of quenches to 
both isomers throughout much of the solid-like region 
until the phase-change region, when the higher-energy 
isomers become accessible. For TV = 57, there are several 
lower-energy isomers below AV m i n ~ 1.0 that were rep- 
resented in the quench distributions, but as with TV = 52 
and 53, there was no corresponding lower-temperature 
heat capacity peak. Both the TV = 57 and 59 quench 
curves show evidence of premelting, consistent with their 
corresponding 5(T) curves in Figure njl 



IV. CONCLUSIONS 

Compared to the irregularity with respect to cluster 
size TV that characterized the properties of the-small clus- 
ters of size TV < 24 that I studied previouslyU the prop- 
erties for the clusters sized 25 < TV < 60 showed more 
regularity and consistency. The lowest-energy isomers 
for these sizes were all dominated by icosahedral based 
geometries, except for TV = 38, whose lowest-energy iso- 
mer had fee geometry (the lowest-energy icosahedral-like 
isomer for N = 38 was only slightly less stable). For 
the icosahedral-like clusters, configurations having anti- 
Mackay overlayers were the lowest-energy isomers for 
N < 31, while configurations having Mackay overlayers 
were the lowest-energy isomers for the remaining sizes up 
to TV = 55, after which the filling of the third shell began 
again with anti-Mackay overlayers. Correspondingly, the 
heat capacity peak parameters formed two overlapping 
sequences as functions of cluster size, with those clus- 
ters having anti-Mackay overlayer lowest-energy isomers 
showing a gradual evolution in their heat capacity peaks 
to higher temperature and smaller size, while those clus- 
ters having Mackay overlayer lowest-energy isomers had 
small, lower-temperature peaks that generally shifted to 
higher temperature and grew in size as TV increased. 

For cluster sizes 25 < TV < 35, the heat capacity 
peak temperatures did not correspond well to the clus- 
ter melting temperature ranges inferred from the tem- 
perature dependence of other cluster properties, such as 
the rms bond length fluctuations and the distributions 
of quenched isomers. For each of the cluster sizes from 
TV = 36 to 60, though, the melting temperature ranges 
obtained from the three measures were in qualitatively 
good agreement. The general trend in the cluster melt- 
ing temperatures as a function of the cluster size over the 
range studied was an initial decrease, reaching a mini- 
mum near TV = 36, then an increase until a maximum at 
TV = 55, where the second Mackay shell was completed, 
followed by another decrease as atoms were added to the 
third Mackay shell. Superimposed on this general trend 
were local maxima corresponding to the magic numbers 



TV = 26, 29, 32, 36, 39, 43, 46, 49 and 55. The heat ca- 
pacity peak heights had a magic number series that was 
a subset of this sequence, with local maxima occuring at 
TV = 36, 39, 43, 46, 49 and 55; for the other magic number 
sizes where the heat capacity peaks did not show magic 
number behavior (TV = 26, 29, and 32), the heat capac- 
ity peaks were located well away from the solid-liquid 
phase-change regions. 

The magic number sequence for the cluster melting 
temperatures coincided very well with the magic num- 
ber sequences obtained for structural properties such as 
the binding energy difference AEb(N), the second finite 
difference of the energy A%E(N), and the gap between 
the lowest-energy and next lowest-energy isomers. The 
only differences were TV = 38 being a magic number for 
the AE b (N) sequence instead of TV = 39, and TV = 58, 
which was a magic number for all three structural prop- 
erties, but was not for either the heat capacity or the rms 
bond length fluctuation sequences. This similarity in the 
magic number sequences is consistent with the widely 
held view that the thermodynamic behavior of clusters 
in the phase-change region is dominated by their poten- 
tial energy surfaces. 

The parallel tempering method was found to be es- 
pecially useful for overcoming quasiergodicity and accu- 
rately calculating the heat capacities of these clusters, 
particularly for the range TV = 30 to 38, where the 
Metropolis method was hard pressed to provide even 
qualitatively correct results. Even the J-walking method, 
which had worked very well for smaller cluster sizes, was 
unable to determine heat capacities reproducibly for clus- 
ters in this range. Given also the relatively small compu- 
tational overhead for the method compared to Metropolis 
Monte Carlo, parallel tempering currently stands as the 
method of choice for cluster Monte Carlo simulations. 
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TABLE I. Cluster heat capacity peak parameters for 25 < N < 60. These values were obtained by smoothing and inter- 
polating the parallel tempering data shown in Figs, o, M and H, which were based on the combined results of 10 walks, each 
consisting of 10 6 passes of data accumulation per temperature. The uncertainties were estimated as the standard errors of the 
set of peak parameters that were obtained from the interpolated heat capacity curves for each of the 10 walks. 
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erature peak 


Higher temperature peak 
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T" 


(C' v )/N 


J" 


(Cv)/N 
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25 






0.3134±0.0016 


6.425±0.011 


43 


26 






0.3116±0.0015 


6.380±0.012 


44 


27 


0.0183±0.0004 


2.965±0.001 


0.3157±0.0015 


6.313±0.011 


45 


28 






0.3242±0.0016 


6.210±0.013 


46 


29 






0.3299±0.0025 


6.092±0.011 


47 
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0.0554±0.0010 


3.059±0.002 


0.3337±0.0020 


5.995±0.014 


48 
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0.0279±0.0019 


3.443±0.079 


0.3396±0.0018 


5.888±0.013 


49 
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0.0603±0.0008 


4.115±0.035 


0.3484±0.0023 


5.814±0.009 


50 
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4.336±0.063 


0.3511±0.0033 


5.737±0.010 
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4.202±0.024 


0.3641±0.0021 


5.656±0.005 
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0.1070±0.0014 


4.759±0.031 


0.3771±0.0064 


5.600±0.010 
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0.4097±0.0074 
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11.702±0.033 

10.934±0.053 

10.340±0.073 
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FIG. 1. Effects of the constraining radius R c on cluster reduced heat capacities per particle for simulations of Ar3o (at left) 
and Ar32 (at right). The filled circles represent parallel tempering results for R c = 3.0 a, while the squares represent parallel 
tempering results for R c = 3.5 a; these were obtained from walks of 10 7 total passes per temperature point. The error bars 
are two standard deviations on each side. The dotted line and the diamonds represent Metropolis results for R c = 4.0 a, 
which were obtained from walks of 10 7 and 10 8 total passes per temperature point, respectively. The curves with different 
constraining radii are in very good agreement on the low-temperature side of the peak, but show substantial deviations on the 
high-temperature side of the peak, effecting the peak height and location slightly. The differences between the R c = 3.5 a and 
R c = 4.0 a peak height and location were minor in these two cases, and insignificant for most of the other cluster sizes whose 
parallel tempering results were obtained using R c — 3.5 a. The smaller lower-temperature peaks seen in the parallel tempering 
curves are absent in the Metropolis curves, indicating the Metropolis simulations were subject to substantial quasiergodicity. 

FIG. 2. Cluster local potential energy minima relative to the global minimum, in reduced units. Magic number behavior can 
be seen for N = 26, 29, 32, 36, 39, 43, 46, 49 and 55, which have large gaps between their global minimum energy and their next 
higher energy configurations, compared to their immediate neighbors. The potential energies were obtained from J-walking 
and parallel tempering quench studies. 

FIG. 3. The binding energy difference as a function of cluster aggregate size (top plot), the second finite difference of 
the energy as a function of cluster size (middle plot), and the energy gap between the lowest-energy isomer and the next 
lowest-energy isomer, as a function of cluster size (bottom plot), all in reduced units. Except for TV = 38 and 39, the three 
measures show the same peak values. The binding energy differences and the second finite differences for 5 < TV < 24 have 
been included for comparison, using data from Ref. [ill. 

FIG. 4. The second finite difference of the energy A2E as a function of temperature T and size TV, for Ar clusters. The 
structure of the zero temperature curve A2-E(TV) (also shown in Fig. H) can be seen to extend well into the higher-temperature 
regions for the magic number sizes (seen as ridges in the surface), which indicates the persistence of their solid- like lowest-energy 
isomers over relatively large temperature ranges. The flat regions of the surface at higher temperature correspond to the cluster 
liquid-like regions, and so the phase-change regions correspond to the temperatures where the ridges abate. 

FIG. 5. Cluster potential energies per particle as functions of temperature, obtained from parallel tempering simulations for 
clusters ranging from TV — 25 (top curve) to N — 60 (bottom curve). Reduced units have been used; the absolute temperature 
scale at the top is for Ar. Magic number sizes for TV = 26, 29, 32, 36, 39, 43, 46, 49 and 55 are indicated by the dotted 
curves, and correspond to those clusters having relatively large energy differences from their preceding cluster sizes at low 
temperatures. A rough indication of the phase-change region from solid-like to liquid-like can be ascertained by the change 
from the irregular spacing between the curves seen at lower temperatures to the regular spacing seen at higher temperatures. 

FIG. 6. Reduced heat capacities per particle as functions of temperature for Ar clusters of size TV — 25 to 36. The open 
circles and thicker solid line represent the results obtained from parallel tempering runs of 10 7 passes per temperature point. 
The dotted line represents the results obtained from Metropolis runs of 10 7 passes per temperature point, while the thin solid 
line represents the results obtained from similar Metropolis runs of 10 8 passes per temperature point. Beginning with Ar3o, 
a smaller, lower-temperature peak can be seen evolving in the parallel tempering curves and moving to higher temperatures, 
while the higher-temperature peak gradually disappears. Quasiergodicity is evident in the Metropolis results of these cluster 
sizes by the absent and malformed lower temperature peaks. As in Fig. |lj, the discrepancies between the parallel tempering 
and Metropolis curves for higher temperatures are due to the parallel tempering simulations having been run with a smaller 
constraining radius of 3.5 a, instead of the value of R c — 4.0 a used in the Metropolis simulations. The large ticks on the 
temperature axes indicate the temperatures where the Metropolis rms bond length fluctuations rose sharply for walks of 10 7 
passes per temperature (higher temperature) and 10 8 passes. 

FIG. 7. Same as Fig. H, but for clusters ranging in size from TV = 37 to 48. Note the different heat capacity scale. 

FIG. 8. Same as Fig. |6j, but for clusters ranging in size from TV — 49 to 60. For these sizes, the constraining radius was set 
to 4.0 a for both the parallel tempering and Metropolis simulations. Again, note the different heat capacity scale. 
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FIG. 9. Magic number behavior in cluster heat capacities. The upper plot shows the reduced heat capacity peak values as 
a function of the cluster size, while the lower plot shows the corresponding peak reduced temperatures as a function of the 
cluster size; the right-hand axis is representative of Ar. The peak parameters were obtained from interpolations of the parallel 
tempering data; parameter values for TV < 25 have been included from Ref. hi] for added context. Two general sequences have 
been identified as a consequence of the clusters in the range 30 < TV < 37 having two prominent peaks. The circles indicate 
the heat capacity parameters for the broader, higher temperature peaks, which have evolved from the solid-liquid phase-change 
region peaks of the smaller cluster sizes, while the squares represent the narrower, lower temperature peaks, which evolve into 
the solid-liquid phase-change region peaks of the larger cluster sizes. The triangles represent the parameters for the very small, 
broad peaks seen for sizes TV = 27 and 58. The magic number labels indicate the cluster sizes whose peaks are higher than 
those of their immediate neighbors. 

FIG. 10. Root mean square bond length fluctuations as functions of the temperature for Ar clusters of size TV = 25 to 36. 
The solid circles represent data obtained from 10 Metropolis walks, each consisting of 10 6 passes per temperature point, while 
the open circles represent data obtained from 10 similar Metropolis walks of 10 7 passes in length. The error bars are single 
standard deviations for the 10 walks. The large ticks on the temperature axes indicate the corresponding cluster heat capacity 
peak temperatures. 

FIG. 11. Same as Fig. lid, but for clusters ranging in size from TV — 37 to 48. 

FIG. 12. Same as Fig. lid, but for clusters ranging in size from TV = 49 to 60. 

FIG. 13. Magic number behavior in cluster rms bond length fluctuations. The curves plot the reduced temperatures associ- 
ated with the Lindemann-like threshold 5 = 0.2, as a function of cluster size. The squares represent the values corresponding 
to S(T) curves whose data were obtained from Metropolis walks of 10 6 passes per temperature point, while the circles represent 
values obtained from similar Metropolis simulations of 10 7 passes. The right-hand axis is representative of Ar. The magic 
numbers are a subset of those obtained in Figure H for the second finite differences in the energy. 

FIG. 14. Quenched isomer distribution curves for Ar clusters as functions of temperature for the magic number sizes 
TV = 36, 39, 43, 46, 49, and 55. The data were obtained from parallel tempering simulations of 10 6 total passes per tem- 
perature. The filled circles represent the fraction of isomers that quenched to the lowest-energy isomer, the open circles the 
second lowest-energy isomer, and the squares the third lowest-energy isomer, while the line represents the sum of the remaining 
isomers. The dotted vertical lines in each plot indicate the corresponding heat capacity peak temperatures. 

FIG. 15. Same as Figure Il4l but for the non-magic number sizes TV = 25, 27 and 31 to 34. 

FIG. 16. Same as Figure [lj, but for the non-magic number sizes TV = 37, 38, 41, 42, 45 and 48. 

FIG. 17. Same as Figure Il4, but for the non-magic number sizes TV = 50, 52, 53, 56, 57 and 59. 
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